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We calculate the pair correlation function of an interacting Bose gas in a harmonic trap 
directly via Path Integral Quantum Monte Carlo simulation for various temperatures and 
compare the numerical result with simple approximative treatments. Around the critical 
temperature of Bose-Einstein condensation, a description based on the Hartree-Fock 
approximation is found to be accurate. At low temperatures the Hartree-Fock approach 
fails and we use a local density approximation based on the Bogoliubov description 
for a homogeneous gas. This approximation agrees with the simulation results at low 
temperatures, where the contribution of the phonon-like modes affects the long range 
behavior of the correlation function. Further we discuss the relation between the pair 
correlation and quantities measured in recent experiments. 
PACS numbers: 03.75.Fi, 02.70.Lq, 05.30.Jp 



I. INTRODUCTION 

One of the appealing features of the experimental achievement of Bose-Einstein condensation in dilute 
vapors , is the demonstration of first order coherence of matter waves . The interference pattern of 
this experiment agrees with the theoretical calculation || , which reveals that the underlying theoretical 
concept of off-diagonal long range order due to a macroscopically occupied quantum state is justified ||. 
Additional experiments have explored certain aspects of second and third order coherence of a trapped 
Bose gas @-^]. Here we study the density-density correlation function which is related to second order 
coherence. With the knowledge of this pair correlation function, the total interaction energy can be 
calculated. In the release energy of the atoms was measured after switching off the magnetic trap. In 
the Thomas Fermi regime at zero temperature the initial kinetic energy can be neglected and the release 
energy is dominated by the interaction energy. By comparison with the usual mean field interaction 
energy using a contact potential, it was concluded that the release energy is mainly proportional to the 
pair correlation function at vanishing relative distance. Strictly speaking this statement cannot be correct 
as for interactions with a repulsive hard core the pair correlation function must vanish at zero distance. 
To give a precise meaning to this statement one needs to access the whole correlation function. 

In this paper we consider in detail the spatial structure of the correlation function of an interacting 
trapped Bose gas. The Fourier transform of this function is directly related to the static structure factor 
which can be probed by off-resonant light scattering. The tendency of bosonic atoms to cluster together 
causes atom-bunching for an ideal gas above the condensation temperature, for the atoms separated by 
less than the thermal de-Broglie wavelengt h ]To|]. For the condensate atoms, this bunching vanishes, since 
they all occupy the same quantum state |1 jjfljfl. However, for a gas with strong repulsive interatomic 
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interaction, it is impossible to find two atoms at exactly the same place, and hence the pair correlation 
function must vanish at very short distances. This mutual repulsion can significantly reduce the amount 
of bosonic bunching at temperatures around the transition temperature |f3[ ]. At much lower temperature, 
the presence of the condensate changes the excitation spectrum as compared to the noninteracting case. 
It is known that in a homogeneous Bose gas the modes of the phonons give rise to a modification of the 
long range behavior of the correlation function |l4j . 

Using path integral quantum Monte Carlo simulations all equilibrium properties of Bose gases can be 
directly computed without any essential approximation [fl5| . ft has been shown that this calculation can 
be performed directly for the particle numbers and temperatures of experimental interest |]l6| . Here, we 
use this approach to calculate the pair correlation function for various temperatures and compare our 
results with simple approximate treatments. 

Near the critical temperature our data are quantitatively well explained by an improved semiclassical 
Hartree-Fock theory, where the full short range behavior is taken into account. At low temperature this 
single-particle approximation fails since the low lying energy modes become important and they are not 
correctly described by the Hartree-Fock treatment. In the Bogoliubov approach these modes are phonon- 
like and change the behavior of the correlation function. Adapting the homogeneous Bogoliubov solution 
locally to the inhomogeneous trap case we find an excellent agreement with the Monte Carlo simulation 
results at low temperature. 



II. HAMILTONIAN OF THE PROBLEM 

The Hamiltonian of N interacting particles in an isotropic harmonic trap with frequency lu is given by 

N r 2 i i i 
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where V is the interatomic potential, which depends only on the relative distance r%j — |r,: — fj \ between 
two particles. This potential in the experiments with alkali atoms has many bound states, so that the 
Bose-condensed gases are metastable systems rather than systems at thermal equilibrium. To circumvent 
this theoretical difficulty, we have to replace the true interaction potential by a model potential with no 
bound states. 

This model potential is chosen in a way that it has the same low energy binary scattering properties as 
the true interaction potential. In the considered experiments, the s-wave contribution strongly dominates 
in a partial wave expansion of the binary scattering problem, so that it is sufficient that the model 
potential have the same s-wave scattering length a as the true potential. For simplicity we take in 
the quantum Monte Carlo calculations a pure hard-core potential with diameter a. In the analytical 
approximations of this paper, we have taken, as commonly done in the literature, the pseudo-potential 
described in fji|] , which is a regularized form of the contact potential, gS(ri — rj) g^-(ri2-), with a 
coupling constant 



4irh 



2, 



111 



(2) 



III. PATH INTEGRAL QUANTUM MONTE CARLO APPROACH 
A. Reminder of the Method 

The partition function Z of the system with inverse temperature /? = (ksT^ 1 is given as the trace 
over the (unnormalized) density matrix g: 

£>(/?) = e^ H (3) 

over all symmetrized states. Both satisfy the usual convolution equation which we can write in the 
position representation: 
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Z = jnT.I d 3N Rg(R,R p ,(3) (4) 
• P J 

= M^/ d3NR J d3NR i-J d 3N R M g(R,R2,r)...g(R M ,R p ,T). (5) 

Here r = /3/M, where M is an arbitrary integer, R is the 3N-dimensional vector of the particle coordinates 
R = (ri,r*2, ...,fjsr), P is a permutation of the AT labels of the atoms and R p denotes the vector with 
permuted labels: R p = (fpfi)) rp(2) > ■■■) ^p(N))- Since only density matrices at higher temperature (r <C /3) 
are involved, high temperature approximations of the A^-body density matrix can be used. 

The simplest approximation is the primitive approximation corresponding to exp[r(j4 + B)] ~ 
exp[ri?/2] exp[rA] exp[r_B/2], which neglects the commutator of the operators A and B. It corresponds 
to a discrete approximation of the Feynman-Kac path integral and gives the correct result in the limit 
M — ► oo [T7|jr^ ]. This can be seen by using the Trotter formula for the exponentials of a sum of two 
noncommuting operators 

f e TA/n e TB/ n y (g) 



e r(A+B) = Um 
n — >oo 



The discretisized path integral for the Af-particle density matrix at inverse temperature r can therefore 
be written in the primitive approximation with symmetric splitting as 



N 

T 



g(R,R , ,T)~Y[g 1 (r k ,f k ',T)Y[e X p[~-(V(r lJ ) + V(f lJ '))\ , (7) 

fc=l i<j 

where gi(fk, t%',t) is the density matrix of nonintcracting particles in the harmonic trap and = r^ — r}, 
r'ij =r' i — r'j. However, this approximation leads to slow convergence since the potential energy in the 
argument of the exponentials are not slowly varying compared to the density matrix of one particle in 
the external potential, Qx(fi, Ti' ,r). This has the consequence that eq.(0) is not a smooth function in the 
region where two particles are in contact, as it should. In order to get such a smooth function we use the 
fact that the potential energy part of eq.(^) can also be written as: 
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(8) 



where the brackets correspond an average over an arbitrary distribution of f*y (t) , starting from and 
ending at r\j ', which reproduces the correct high temperature limit of the primitive approximation. It 
is convenient to take the random walk corresponding to the kinetic energy as weight function so that 172 
is the solution of the binary scattering problem in free space: 



^ ~, ,_ gjj exp[-T(p?./m + V(r l3 ))M 3 ) (Q , 
(ry|exp[-T^/TO]|rU 



where pij is the operator of the relative momentum between particles i and j. This leads to the so called 
pair-product approximation |l8| , |l5| ], where the density matrix is approximated as 

JV 

q(R,R',t)~ Yl ei(*n,»V,T) Y[92{f%j\fij',r). (10) 

n— 1 i<j 

This approximation has the advantage to include exactly all binary collisions of atoms in free space, only 
three and more atoms in close proximity will lead to an error; convergency with respect to M — > 00 is 
reached much faster. In the simulation the two-particle correlation function §2 is equal to one for non- 
interacting particles and plays the role of a binary correction term in presence of two-body interactions. 

As in |1(|] we take A^ = 10,000 particles with a hard-core radius of a — 0.0043(?i/mct;) 1 / 2 . The 
transition temperature of the noninteracting Bose-gas is hsT® = 20.26 hui or /3° ~ 0.05(?iaj) -1 and a 
value of r = 0.01(?iw) _1 was found sufficient. In the low temperature regime (kgT -C fi 2 /ma 2 ) the 
most important contribution to (72 for hard spheres is the s-wave contribution, which can be calculated 
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analytically |l9[; for non vanishing relative angular momenta (I > 0) we neglect the effect of the potential 
outside of the hard core. In this way we can obtain an explicit formula for g 2 , 



g 2 (r,r";T) = 1 + 



h 2 [3 



1 



exp 



(r + r') 



>\2 



exp 



(r + r' - 2a) 2 
Ah 2 (3/m 



-(rr f -\-r-f f )rn/2H 2 (3 



(ii) 



4?T/3/m. 

for r and r' outside of the hard core diameter (|r| > a and \ f'\ > a), otherwise g 2 = 0. 
The density-density correlation function can be easily calculated as 

As the atoms are in a trap rather than in free space, this quantity is not a function of the relative 



(12) 



coordinates 



of the two particles only. Imagine however that this pair distribution function 



be probed experimentally by scattering of light by the atomic gas, where we assume a large beam 
waist compared to the atomic sample. As the Doppler effect due to the atomic motion is negligible, 
the scattering cross section depends only on the spatial distribution of the atoms. Furthermore, for 
a weak light field very far detuned from the atomic transitions, the scattering cross section can be 
calculated in the Born approximation; it then depends only on the distribution function of the relative 
coordinates f ' — r" between pairs of atoms. We therefore take the trace over the center-of-mass position 
R = (f'+r")/2: 

V i2) (r,P)= *_ X) J dW 2) (i? + r72;i?-r/2,/3), (13) 

where we have divided by the number of pairs of atoms to normalize (p^ to unity. Note that the result 
depends only on the modulus r of r as the trapping potential is isotropic. 



B. Results of the Simulation 

In fig.l we show j3) for various temperatures below T®, obtained by the simulation of the 

interacting bosons in the harmonic trap, where the critical temperature T c is reduced compared to the 
ideal gas |2^,|l^,|l[ . All pair correlation functions are zero in the region of the hard-core radius as they 
should. At larger length scales the r dependence of the result is also simple to understand qualitatively, 
as we discuss now. 

Consider first the case T > T C1 where no condensate is present. As the typical interaction energy n(r)g 
(n(r) being the total one-particle density at r) is much smaller than fceT, we expect to recover results 
close to the ideal Bose gas. The size of the thermal cloud (ksT /muj) 1 / 2 determines the spatial extent of 
^ 2 \r); the bosonic statistics leads to a spatial bunching of the particles with a length scale given by the 
thermal de Broglie wavelength 



' 2nTi 2 
inksT 



(14) 



The Bose enhancement of the pair distribution function is maximal and equal to a factor of 2 for particles 
at the same location (r = 0). This effect is preserved by the integration over the center of mass variable 
and manifests itself through a bump on ip^ 2 \r) in fig.l. Due to the influence of interactions the bump is 
suppressed at small distances and the factor of 2 is not completely obtained. 

For T < T c a significant fraction of the particles accumulate in the condensate. As the size of the 
condensate is smaller than that of the thermal cloud, the contribution to of the condensed particles 
has a smaller spatial extension, giving rise to wings with two spatial components in ip( 2 \ as seen in fig.l. 
Apart from this geometrical effect the building up of a condensate also affects the bosonic correlations at 
the scale of At: The bosonic bunching at this scale no longer exists for particles in the condensate. This 
property, referred to as a second order coherence property of the condensate @,§.|l3|, is well understood 
in the limiting case T = 0; neglecting corrections due to interactions, all the particles are in the same 
quantum state \tpo) 80 that e.g. the 2-body density matrix factorizes in a product of one-particle pure 
state density matrices. This reveals the absence of spatial correlations between the condensed particles. 
This explains why in fig.l the relative height of the exchange bump with respect to the total height is 
reduced when T is lowered, that is when the number of non-condensed particles is decreased. 



4 



IV. COMPARISON WITH SIMPLE APPROXIMATE TREATMENTS 



At this stage a quantitative comparison of the Quantum Monte Carlo results with well known approx- 
imations can be made. 



A. In presence of a significant thermal cloud: Hartree-Fock approximation 

As shown in |2l|| in detail, at temperatures sufficiently away from the critical temperature, the Hartree- 
Fock approximation pp|| gives a very good description of the thermodynamic one-particle properties. 

To derive the Hartree-Fock Hamiltonian we start from the second quantized form of the Hamiltonian 
with contact potential 



H= d 



j3 ,- 



&(r)(H + |# + (^* + (f)#(f)*(f) ( 15 ) 

where Hq is the single particle part of the Hamiltonian. Due to the presence of the condensate we split 
the field operator $ in a classical part i/jq, corresponding to the macroscopically occupied ground state 
and the part of the thermal atoms tj) with vanishing expectation value (ijj) = 0: 

#(r) ~ Vo(f) + V>(r). (16) 

After this separation we make a "quadratization" of the Hamiltonian by replacing the interaction term 
by a sum over all binary contractions of the field operator, keeping one or two operators uncontracted, 
e.g. 

V^t^i ~ - 2(?// t ?/>)(V> t V>). (17) 

This is done in such a way that the mean value of the right hand side agrees with the mean value of 
the left hand side in the spirit of Wick's theorem. In the Hartree-Fock approximation we neglect the 
anomalous operators, such as ■0 T, T , and their averages, and we end up with a Hamiltonian which is 
quadratic in ipo and V>, but also linear in i[j and ijy . Now we choose tpo such that these linear terms vanish 
in order to force = 0. This gives the Gross-Pitaevskii equation for the condensate p2| 

H 2 \7 2 1 1 

— \- -muj 2 r 2 + g[n (r) + 2n T (r,f)] > i/Jo(r) = /#oM (18) 

2m 2 J 

where no(r) = |-0o(r)| 2 corresponds to the condensate density with Nq particles and nr(r, 'f) — 
{r)ip{r)) is the density of the thermal cloud. 
Up to a constant term we are left with the Hamiltonian for the thermal atoms 

H = J d 3 rft(r){H + 2gn{r)- ^)^(r) (19) 

where nir) — n (r) + u-r(f, r) denotes the total density and depends only on the modulus of r. To work 
out the density-density correlation function, we formulate ([l2]) in second quantization: 

g (2) (r;r',f3) = (r)¥ (r')V(r ')*(r)), (20) 

we use the splitting Jl6|), together with Wick's theorem and get 

g^ F {f;r',P) = Mr)Mr)Mr')Mr') 

+ipo(r)ip (r)n T (r', r') + ip (r')ip (r')n T (r, r) + 2ij^(r)^ (r')n T (r, r') 

+n T {r 1 f)n T {r' , r') + n T {r, r")n T {r, r"). (21) 

Here we have chosen the condensate wave function to be real and 

n T (r ) r') = (ft(f)iP(f')) (22) 



5 



corresponds to the nondiagonal elements of the thermal one body density matrix. Since the Hamiltonian 
( |l9| ) of the thermal atoms is quadratic in ip, this density matrix is given by 



riT(f, f ') = (r\ 



exp \ft(H + 2gn(r) - fi)] - 1 1 



(23) 



In the semiclassical approximation (fc^T 3> tku) we can calculate explicitly these matrix elements by 
using the Trotter break-up, which neglects the commutator of r and p: 



n T (r,r r ) = ^^e-Wfe+W^+a^W-/*)^/) 
1=1 



(24) 



We finally get 



1 = 1 



(25) 



1 °° 1 



p/2 



exp 



f/|2 



IX'j-i 



1/3 (muj 2 (r 2 + r' 2 )/4 + g{n(r) + n(r')) - fj) 



(26) 



For the diagonal elements the summation gives immediatly the Bose function ff 3 / 2 (g) = Yli^i z l /l 3 ^ 2 - 
For a given number of particles N, eq.(|l^) and the diagonal elements r = f'of eq.(p6|) have to be solved 
self consistently to get the condensate density no(r) and the thermal cloud riT{r,r)- With this solution 
we can work out the nondiagonal matrix elements of the density operator which give rise to the exchange 
contribution of the density-density correlation (21), and the correlation function can be written as a sum 
over the direct and the exchange contribution 



q {2) , ir;r'0). 

t-exenange^ ' 



(27) 



Up to now the short range correlations due to the hard core repulsion have not been taken into account, 
but we can improve the Hartree-Fock scheme further to include the fact that it is impossible to find two 
atoms at the same location: We assume that the particle at r interacts with the full Hamiltonian with the 
particle at r ' but only with the mean-field of all others (over which we integrated to derive the reduced 
density matrix). This gives in first approximation: 



f 2) F (r- r ', f}) = g {2 l ect (r; f, f3)g 2 (r~ r>; f- r ', 0) + Q™ change (r; f ', P)g 2 (r- 



r ; r 



f,/3) 



Q [ HF(?-,r' ,P)g2{r- f';f- r',/3) 



(28) 
(29) 



where the two particle correlation function g 2 is the solution of the binary scattering problem, eq.(|ll[). 
Further we used the fact that g 2 — 1 for particle distances of the order of At and larger. In principle 
one should integrate over the second particle to get a new one-particle density matrix and find a self- 
consistent solution of the Hamiltonian. But since the range of g2 is of the order of the thermal wavelength, 
it will only slightly affect the density, so we neglect this iteration procedure. Using the solution of the 
coupled Hartree-Fock equations to calculate (p9|), and integrating over the center-mass-coordinate, we 
get ip HF (r, @). As shown in fig.l, this gives a surprisingly good description of the correlation function at 
high and intermediate temperatures. 
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FIG. 1. Pair correlation function (p^ (r, /3) vs r in units of the harmonic oscillator length (h/mui) 1 / 2 from the 
Monte-Carlo and the Hartree-Fock calculations for (3 = 0.05(ftu) _1 , (3 = 0.06(ftu) _1 , and (3 = 0.07(Rw) _1 (from 
the bottom to the top). The corresponding condensate fractions N /N are 0.0 (T : /T° ~ 1.0), 0.22 (T/T c ° ~ 0.8), 
and 0.45 (T/T® — 0.7). T,? is the Bose-Einstein condensation temperature for the ideal gas. For clarity we removed 
the part of short r for the upper curves. 



B. The quasi-pure condensate: Bogoliubov approach 



The Hartree-Fock description must fail near zero temperature: Since the anomalous operators i/j'ij)' 
and iptp have been neglected, it describes not well the low energy excitations of the systems. It is known 
that the zero temperature behavior can be well described by the Bogoliubov approximation j23|. In 
this paper it is not our purpose to calculate the correlation function using the complete Bogoliubov 
approach in the inhomogeneous trap potential. This could be performed using approaches developed 
in |^,^5|. Here we use the homogeneous description of the Bogoliubov approximation and adapt it to 
the inhomogeneous trap case with a local density approximation. This approach already includes the 
essential features which the Hartree-Fock description neglects at low temperatures. 

We start with the description of the homogeneous system with quantization volume V and uniform 
density n = N/V. As in ||(| we split the field operator ^> into a macroscopically populated state <E> and 
a remainder, which accounts for the noncondensed particles: 



#(r) = $(r)a$ + Sif(r). 



(30) 



In the thermodynamic limit N — > oo, V — > oo, keeping N/V — n and Ng — const, the typical matrix 
elements of 5^ at low temperatures are yJ~N times smaller than a$ . Hence we can neglect terms cubic 
and quartic in <5^, when we insert (^0|) in the expression of the density-density correlation function (pp[). 
Since the condensate density is given by the total density minus the density of the excited atoms, we have 
to express the operator of the condensate density in the same order of approximation for consistency, 



1 



W = N - 



Finally we use the mode decomposition of the homogeneous system 



fc#0 



■hie 

k 



—ik-r„ 



(31) 



(32) 



where bj: annihilates a quasiparticle with momentum k. The components itfc and satisfy the following 
equations: 
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2m ' 9n 

-gn 




= E k 



Uk 
Vk 



together with the normalization: 



\Vk? = I- 



(33) 



(34) 



At low temperatures the quasiparticles have negligible interactions and we can use Wick's theorem to 
get the following expression for the correlation function 



QBG(i 



',(3) =n 2 + 2n 



d 3 fce *fe-(r'-r") 



iu\ 



2ufcUfc)(Sl6g / 



- u k v k 



+ 0(Vn) (35) 



where we used <£>(r) = V 1 / 2 . The quasiparticles obey Bose statistics, so that the mean number of 
quasiparticles with momentum k and energy E k is given by 



1 



o pE k 



(36) 



We see from eq. (|3g) that in the homogeneous system the density-density correlation function depends 
only on the relative distance r — \r ' — r"\. The derivation of the properties of the pair correlation function 
is given in the appendix. At T = the pair correlation function has the following behavior |]l4| , ^7f 



9„ 



t(r,T = Q) 



1 



(37) 



(r » 



where £ = (Snna)" 1 / 2 is the healing length of the condensate. For finite but small temperatures this 
structure is only slightly changed (see appendix). The modification of the low energy spectrum due to 
the Bogoliubov approach is responsible for the long range part of the correlation function. 

Apart from the edge of the condensate, the total density n(r) for low temperature in the trapped 
system varies rather slowly compared to the healing length £ for the considered parameters. So it is 
possible to adapt the result of the homogeneous system to the inhomogeneous trap case. For a given 
density n(r) we get with a local density approximation for the pair correlation function instead of eq.(|3T 



^(r,/3)c f d 3 R {n(\R + f/2\)n(\R-f/2\) 



-2n(R) / d 3 ke 



Ak-r 



(u 2 k (R) + v 2 k (R) + 2u k (R)v k (R)) (< 



v 2 k (R) + u k (R)v k (R) 



where Uk(R), v k {R), and E k (R) are solutions of eq.(|33J) for the given density n(R). 

As shown in fig. 2 this gives an excellent agreement with the Quantum Monte Carlo results at low 
temperature. We have checked that at this temperature the difference with the Bogoliubov solution at 
T = is almost negligible. The good agreement with the simulation reflects that the long range behavior 
of the pair correlation function in this approximation is correctly described by eq. ( |37| ) . We note that in an 
intermediate temperature regime, which is not shown, both approaches, the Hartree-Fock and the local 
density Bogoliubov calculation, do not reproduce the simulation results quantitatively: The maximum 
local error is about 5%. 



(38) 
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FIG. 2. Pair correlation function tp^ (r, /3) vs r in units of the harmonic oscillator length (h/mui) 1 / 2 from the 
Monte Carlo, the Bogoliubov and the Hartree-Fock calculations for (3 — 0.20(?itj) _1 with a condensate fraction 
Nq/N ~ 0.95 (T/T c ° ~ 0.25). The healing length is roughly £ ~ 0.3 in this units 



V. CONNECTION TO THE INTERACTION ENERGY 



The knowledge of the pair correlation function permits us to calculate the total energy of the trapped 
atoms E t ot- 



E,, 



1 



Tr{p} 



[Tv{H Q p} + Tv{H lP }} 



(39) 



One has to pay attention that the regularized form of the contact potential, V = g§{r)-§p (r-), acts 
on the off-diagonal elements r±2 and r' 12 of the density-density correlation function Lp(r' 12 ,r 12, (3) = 
(^1 'j ?2 T2) in the space of relative coordinates r±2 and 



12- 



As the 2-body density matrix 



l P{ r i2i r i2j 0) diverges as (1 — a/r' 12 )(l — ajryi) we actually get the simple form: 



1 



Tr{p} 



Tr{£f/p} 



N(N - l)g 



df 



6(r) d 
r dr 



(rV (2) (^/?)) 



(40) 



This form involves only the diagonal elements of the correlation function ip( 2 >(r,/3). Both the improved 
Hartree-Fock solution and the Bogoliubov solution behave for small distances (r <§; £) like 



cpW ( r ~ 0, /3) ~ (1 - a/r) V (2) (0, 0) 



(41) 



where ^ (2) (0,/3) can be obtained graphically by extrapolating the pair correlation function to zero, 
neglecting the short range behavior (r < £); numerically it can be obtained from the Hartree-Fock 
calculation of (^Tj) (see |l3| for analysis of the temperature dependence of <^ 2 ^(0, 0) ). This behavior of 



the correlation functions shows that eq. ( fj 
with the mean interaction energy (Hi): 

(Hi) =s g 



gives a finite contribution linear in a, which we can identify 



N(N-l) ^ (2) 



^(0,/3). 



(42) 



In order g 2 , eq.(^0|) contains a diverging part, We note without proof that this divergency is compensated 
within the Bogoliubov theory by a divergent part of the kinetic energy, so that the mean total energy, 
eq.(|3^), is finite. This lacks in the Hartree-Fock calculation, which is, however, limited to linear order of 

g- 
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In the Thomas-Fermi limit the kinetic energy is negligible, and the interaction energy eq.(f42[) domi- 
nates the total energy, which can be measured. This measurement provides some information about the 
correlation function, however, the true correlation function is not accessible. Only the fictive correlation 
function ip^ (0, j3) for vanishing interparticle distances is obtained. 



VI. CONCLUSION 



We numerically calculated the pair correlation function of a trapped interacting Bose gas with a 
Quantum Monte Carlo simulation using parameters typical for recent experiments of Bose-Einstein con- 
densation in dilute atomic gases. At temperatures around the critical point, an improved Hartree-Fock 
approximation was found to be in good quantitative agreement with the Monte Carlo results. The im- 
proved Hartree-Fock calculation presented in this paper takes the short-range behavior of the correlation 
function into account, especially the fact that two particles can never be found at the same location. 
At low temperature we compared our simulation results to a local density approximation based on the 
homogeneous Bogoliubov approach. The phonon spectrum changes the behavior of the pair correlation 
function for distances r of the order of the healing length £. With the knowledge of the pair correlation 
function we calculated the total interaction energy. We showed that the results of recent experiments 
on second order coherence do not measure the true correlation function, which has to vanish for small 
interparticle distances. Only an extrapolated correlation function is determined, where the exact short 
range behavior disappears. 
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VII. APPENDIX 



In this appendix we give the explicit formulas for the pair correlation function in the Bogoliubov 
approach for an homogeneous system and discuss its behavior at short and long distances, since only 
some aspects have been discussed in literature |L4|,^8|. Starting from eq.(|35|), the pair correlation function 



(2) 

bn=const can be be written explicitly as: 



1 16a f 

L- f ,, ns t( r >P) = 77 ! + / dqsin(qR)f{q) 

f [ f r Jo 



(43) 



with R = y/2r/£, (£ 



iirna) 1 / 2 is the definition of the healing length) and 

-i-i / 



/(?) = 



1 



q 



2 + 



- 1 



(44) 



To get the behavior of eq.(|43|) for small distances (r <C £), we can replace /(g) by its behavior for large 
wavevectors, q — > oo 



/(<?) 



1 



q — > oo. 



Using the value of the integral 1 29 



dx ■ 



SIM 7T 



(45) 



(46) 



we get the short range behavior of the pair correlation function |27|: 
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^n— const (T' 0) 



1 

V 



a~ 


1 


a' 


1 - 2- 




1 - - 


r. 


~ V 


r . 



r < £. 



To get the long range behavior (r>(), we integrate several times by part: 

f°° 1 1 1 

jf d g sm(qR)f(q) = -/(0) - ^/ (2) (0) + ^/ (4) (0) 

For the function f(q) and its derivatives at q — we get 



(47) 



(48) 



T = 0: /(0)=0, / (2) (0) = 1 

T^O: /(0)=0, / (2) (0)=0, / (4) (0)=0, 



(49) 



and the long range behavior at zero temperature given in (|3?]) is obtained. For finite temperature it can 
be shown that /(g) is an odd function of q, so that /^"•'(O) = for all n. Due to that the correlation 
function vanishes faster than any power law in 1/R. 

To work out an explicit expression for finite temperatures we use this antisymmetry to extend the 
range of the integral (^3|) to — oo and we can analytically calculate the expression for two limiting cases 
via the residue calculus. For large distances we only have to take the poles qo of f(q) with the smallest 
modulus into account. For Xt/2tt <C £ corresponding to fc^T ^> ng, and r > we get qo = i, so that 



„ 1 A T 
1 + 2^t; exp 



(50) 



Note the + sign in this expression, leading to <l>n=canst > V^> ^ na * we interpret as a bosonic bunching 
effect for thermal atoms. In the opposite limit, \t/2tt ^ £ and r 3> A T /47r 2 £, the pole with the smallest 
imaginary part is given by go = i4-7r 2 £ 2 /A T and we get [ p8| 



^n— const (Ti ft) y 



(2tt) 3 4< 4 
n A T r 



exp 



-An- 



(51) 



[1] M. H. Anderson, J. R. Ensher, M. R. Matthews, C. E. Wieman, and E. A. Cornell, Science 269, 198 (1995). 
[2] K. B. Davis, M.-O. Mewes, M. R. Andrews, N. J. van Druten, D. S. Durfee, D. M. Kurn, and W. Ketterle, 

Phys. Rev. Lett. 75, 3969 (1995). 
[3] C. C. Bradley, C. A. Sackett, J. J. Tolett, and R. G. Hulet, Phys. Rev. Lett. 75, 1687 (1995); C. C. Bradley, 

C. A. Sackett, and R. G. Hulet, Phys. Rev. Lett. 78, 985 (1997). 
[4] M. R. Andrews, C.G. Townsend, H.-J. Miesner,D.S. Durfee, D.M. Kurn, and W. Ketterle, Science 275, 637 

(1997). 

[5] A. Rohrl, M. Naraschewski, A. Schenzle, and H. Wallis, Phys. Rev. Lett. 78, 4143 (1997). 

[6] C.N. Yang, Rev. Mod. Phys. 34, 694 (1962). 

[7] W. Ketterle and H.-J. Miesner, Phys. Rev. A 56, 3291 (1997). 

[8] E.A. Burt, R.W. Christ, C.J. Myatt, M.J. Holland, E.A. Cornell, and C.E. Wieman, Phys. Rev. Lett. 79, 
337 (1997). 

[9] Yu. Kagan, B.V. Svistunov, and G.V. Shlyapnikov, Pisma. Zh. Eksp. Teor. Fiz. 42, 169 (1985) [JETP Lett. 
42, 209 (1985)]. 
[10] L. Van Hove, Phys. Rev. 95, 249 (1954). 
[11] F. London, J. Chem. Phys. 11, 203 (1942). 

[12] F. Brosens, J.T. Devreese, and L. F. Lemmen s, Phys. Rev. E 55 6795 (1997). 



[13] M. Naraschewski and R. J. Glauber, preprint cond-mat/9806362 



[14] K. Huang, Statistical Mechanics, (John Wiley & Sons, MA 1987), chapter 13; T.D. Lee, K. Huang, and C.N. 

Yang, Phys. Rev. 106, 1135 (1957). 
[15] E.L. Pollock and D.M. Ceperley, Phys. Rev. B30, 2555 (1984); B36, 8343 (1987); D.M. Ceperley, Rev. Mod. 

Phys. 67, 1601 (1995). 



11 



[16] W. Krauth, Phys. Rev. Lett. 77, 3695 (1996). 

[17] R.P. Feynman, Statistical Mechanics (Benjamin/ Cummings, Reading, MA, 1972). 
[18] J.A. Barker, J. Chem. Phys. 70, 2914 (1979). 
[19] S.Y. Larsen, J. Chem. Phys. 48, 1701 (1968). 

[20] V.V. Goldman, I.F. Silvera, and A.J. Leggett, Phys. Rev. B 24, 2870 (1981); S. Giorgini, L.P. Pitaevskii, 
and S. Stringari, Phys. Rev. A 54, R4633 (1996). 



[21] M. Holzmann, M. Naraschewski, and W. Krauth, preprint oond-mat/9806201 



[22] E.P. Gross, Nuovo Cimento 20, 454 (1961); L.P. Pitaevskii, Sov. Phys. JETP 13, 451 (1961). 

[23] N.N. Bogoliubov, J. Phys. (Moscow) 11, 23 (1947). 

[24] A.-C. Wu and A. Griffin, Phys. Rev. A 54, 4204 (1996). 

[25] A. Csordas, R. Graham, and P. Szepfalusy, Phys. Rev. A 57, 4669 (1998). 

[26] C. Gardiner, Phys. Rev. A 56 1414 (1997); Y. Castin, and R. Dum, Phys. Rev A57, 3008 (1998). 

[27] For the short range contribution r < ^ we get in our approach (1 — 2a/r)/V instead of (1 — a/r) 2 /V, which 

is the true short range behavior. But the correcting term in order a 2 is of higher order than the calculation. 

See also @. 

[28] E.M. Lifshitz and L.P. Pitaevskii, Statatistical Physics, Part 2 (Pergamon Press, Oxford, 1980), Chapter 9. 
[29] M. Abramowitz and LA. Stegun, Handbook of mathematical functions, (Dover Publications, New York, 1972), 
Chapter 5. 



12 



